function [mtest,dmtest,vf,vr,vt,a]=steadystate(m,kf,kr,S,epsilon,tolerance,t)
% generic code intended to work with any network
% this code is not optimized for paralellization and will run much slower
% than needed. The same structure could be used to generate ODEs instead of
% calculating metabolites concentrations using the full stochiometric
% matrices for forward and reverse reactions
% kf and kr are forward rate constants
% Sf and Sr the substrate stoichiometric matrices
% epsilon is a small constant to prevent division by zero, and can also be 
% used to moderate the amount of scaling by "a"

% S is the full stoichiometric matrix

Sf=-1*S.*(S<0);
Sr=S.*(S>0);

n=size(Sf,2);
p=size(Sf,1);

a=ones(p,1);
vf=ones(n,1);
vr=ones(n,1);
dm=ones(p,1);
%m=ones(p,1);
mtest=ones(p,0);
dmtest=ones(p,0);

count=0;

% This routine is run untill steady state is reached. This is determined by
% reaching a setpoint value of the lowest "a" value of all metabolites.
% Because "a" is inversely related to metabolite concentraion change, one
% could either use max(m) or min(a)

%while max(abs(dm))>tolerance
% if your problem is not reaching steady state, replace the while loop with 
% the for loop below to troubleshoot      
 for test=1:1   
    for i=1:n % cycle through all fluxes
        vf(i)=kf(i);
        vr(i)=kr(i);
            % take the product of all metabolite concentrations raised to
            % their stoichiometries. Because the matrix is sparse, most
            % stoichiometries are 0, resulting mainly in inconsequential 
            % multiplications of 1
            for j=1:p
                vf(i)=vf(i)*m(j)^Sf(j,i);
                vr(i)=vr(i)*m(j)^Sr(j,i);
            end
    end
    % in the GPU solution each core should evaluate a single row of S 
    dm=S*(vf-vr);
    %dm=-Sf*vf+Sr*vr;
    % m=m+dm*t;
    m=m+a.*dm*t;
    m=m.*(m>0)+(m<0)*1e-5;
    % a=(a+epsilon+m./(abs(dm)+epsilon))./2;
    a=epsilon+1./(abs(dm)+epsilon);
    mtest=cat(2,mtest,m);
    dmtest=cat(2,dmtest,dm);
    count=count+1;
end
vt=vf-vr;
% plotting your results will help you to figure out what is going on 
subplot(2,1,1); plot(1:count,mtest);
subplot(2,1,2); plot(1:count,dmtest);
    
                    
        